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Spin-AS Heisenberg quantum antiferromagnets on the Kagome lattice offer, when placed in a 
magnetic field, a fantastic playground to observe exotic phases of matter with (magnetic analogs of) 
superfiuid, charge, bond or nematic orders, or a coexistence of several of the latter. In this context, 
we have obtained the (zero temperature) phase diagrams up to S' = 2 directly in the thermodynamic 
limit thanks to infinite Projected Entangled Pair States (iPEPS), a tensor network numerical tool. 

We find incompressible phases characterized by a magnetization plateau vs held and stabilized 
by spontaneous breaking of point group or lattice translation symmetry(ies). The nature of such 
phases may be semi-classical, as the plateaus at |th, (1 — ^)th and (1 — ^)th of the saturated 
magnetization (the latter followed by a macroscopic magnetization jump), or fully quantum as 
the spin-1 |-plateau exhibiting coexistence of charge and bond orders. Upon restoration of the 
spin rotation U{1) symmetry a hnite compressibility appears, although lattice symmetry breaking 
persists. Eor integer spin values we also identify spin gapped phases at low enough held, such as the 
S = 2 (topologically trivial) spin liquid with no symmetry breaking, neither spin nor lattice. 

PACS numbers: 


Introduction. - The antiferromagnetic quantum Heisen¬ 
berg model on the Kagome lattice (KHAF) is one of the 
most intriguing strongly-correlated systems. Because of 
the geometric frustration of the lattice, the spin-| case 
{S = ^) has proven exceptionally hard to understand: af¬ 
ter more than twenty years of study, there is still no clear 
consensus about the nature of its ground state [1], though 
recent Density Matrix Renormalisation Group [2] and 
Variational Monte Carlo [3] computations suggest the 
emergence of a topological or critical quantum spin liquid 
(SL), respectively. Tensor network numerical tools have 
come up with good variational energies for the KHAF [4] . 
In the presence of a magnetic held the system exhibits 
several incompressible phases detectable as plateaus in 
the longitudinal magnetization [5, 6]. An external mag¬ 
netic held may also give rise to additional plateaus which 
behave like chiral spin liquids in the XY regime [7]. The 
model has also been studied for larger values of the spin, 
for which the strong quantum huctuations present in the 
S = ^ case are weaker. For instance, for spin S = 1 
the ground state (GS) is believed to be a simplex solid 
[8, 9] (yet other alternatives have also been put forward 
[10-12]), and in the presence of a held a rich phase 
diagram was predicted including nematic and superne¬ 
matic phases [9] . It is also suspected that incompressible 
phases for non-zero held persist for larger spin values [13]. 
The KHAF is also of experimental relevance, e.g., miner¬ 
als such as herbersmithite ZnGu 3 (OH) 6 Cl 2 , volborthite 
Gu 3 V 207 ( 0 H )2 • 2 H 2 O and vesignieite BaGu 3 (OH) 6 Cl 2 
can be described by the S = ^ KHAF [14]. Some 
nickelate or vanadate compounds such as KV 3 Ge 209 
and BaNi 3 ( 0 H) 2 (V 04)2 consist also of weakly coupled 
S = 1 KHAF layers [15], and minerals such as chromium- 
jarosite KGr 3 ( 0 H) 6 (S 04)2 have a KHAF structure with 
spin S' = I [13]. 

In fact, the KHAF is a fantastic toy model to study 


the crossover between quantum and classical mechanics. 
This is so because for the smallest-possible quantum spin 
value, i.e., S = |, quantum huctuations due to geometric 
frustration are so strong that the GS at zero held stabi¬ 
lizes, to our best evidence so far, in a quantum spin liquid, 
which is a inherently quantum-mechanical phase. How¬ 
ever, as the spin gets larger, quantum huctuations get 
weaker as the classical limit S' ^ oc is approached. Un¬ 
derstanding this quantum-to-classical cross-over would 
throw some more light on how quantum ehFects can be 
enhanced or suppressed in quantum lattice systems. In 
turn, this may also provide guidance in the search of new 
exotic phases of quantum matter, and novel quantum 
materials. 

Motivated by this, here we analyze the zero- 
temperature phase diagram of the KHAF up to spin 
S = 2 with a magnetic held. For the numerical simu¬ 
lations we use tensor network methods [16] based on in- 
hnite Projected Entangled Pair States (iPEPS) [17, 18], 
which work directly in the thermodynamic limit. In the 
presence of a held we hnd a host of quantum phases 
of matter at zero temperature for dihFerent spin values. 
These phases may correspond to a magnetization plateau 
(incompressible) or not (compressible), and also break 
spontaneously or not U{1) (spin rotation), translation, 
Cq and/or rehexion symmetries, see Table I. The as¬ 
sociated long range orders are the magnetic analogs of 
the electronic charge density (here named “Solid”), bond 
density (here named “Valence Bond Grystal” or VBG), 
nematic and superhuid orders. Incompressible phases, 
corresponding to plateaus, shrink quickly - or disappear 
- when the value S of the spin increases. Moreover, 
we also hnd macroscopic magnetization jumps at large 
helds, in accordance with a description of independent 
magnons.[19] 

Model and methods.- Here we consider the antiferro- 
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FIG. 1: [Color online] Two patterns from the 9 sites unit cell 
calculations which (a) preserves (the original 3 sites unit cell 
is recovered) or (b) breaks translation symmetry (\/3 x \/3 
superstructure). In the figure, the cr^ and <7v reflexion sym¬ 
metry lines are represented by yellow and green dashed lines, 
respectively, ^-rotations around the crossing point of the 
two dashed lines generate a Cq symmetry. We also show 
the phases (c) Solid 1, (d) Solid 2 and (e) VBC-Solid from 
Table I. The circles represent the longitudinal magnetization 
{S^) at every site, and the lines represent the local energy 
terms (Si • Sj). In all cases, red (blue) means positive (neg¬ 
ative), with the size/thickness proportional to the absolute 
value. These three plots correspond, in particular, to the 
three plateaus in the aS = | phase diagram at magnetization 
I (c), I (d) and | (e) in Fig. 2(a). Notice that the phase 
(e) is, conceptually, a VBC with coexisting solid order, rather 
than a pure solid phase. If the valence bonds strongly res¬ 
onate, the VBC order “melts” and one recovers the Solid 2 
phase in (d). 


magnetic quantum Heisenberg model in the presence of 
an external magnetic field h along the 2 ; direction on the 
Kagome lattice (KHAF), 

H = Y,Si-Sj-hY,St. ( 1 ) 

{ij) « 

We will use values of the spin S' up to S' = 2. Note that 
spin-rotation around z {U{ 1 ) symmetry) is preserved. As 
discussed above, this model exhibits a rich variety of ex¬ 
otic quantum behaviours for different spin values. 

Our numerical calculations are done using tensor net¬ 
work methods based on iPEPS [17]. More specifically, 
we run algorithms for imaginary-time evolution in order 
to obtain approximations of the GS of the system in the 
thermodynamic limit. As explained in the supplementary 
material, we use two different algorithms. Algorithm 1 is 
based on mapping the kagome lattice to the square lattice 
while keeping the locality of interactions, and applying 
then the simple update for a square-lattice PEPS. We use 
such method to estimate the ground-state energy at zero 
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TABLE I: Comparison between different held-induced phases 
as characterized by their preserved (/) or broken (x) point 
group Ce, translation group Gt, spin-f/(l) and reflexion cr^ 
and ay symmetries. For the ^T-symmetry breaking, the phase 
is depicted by 9 sites per unit cell, the so-called a/3 x a/3 
phase. The rotation symmetry Cq can be broken either into a 
lower rotation symmetry group C3 or C2, or completely (x). 
Parameter d is the GS degeneracy, in accordance with the 
remaining discrete symmetries. The spin liquid listed here 
is topologically trivial (since d = 1). When the continuous 
U{1) symmetry is broken, a zero-energy (Goldstone) mode is 
expected. 


field for the spin one-half case (see also supplementary 
material). In algorithm 2, we use a simplex representa¬ 
tion of the tensor network [4] and a simple update scheme 
as well [18], and given its efficiency, we use it for the 
finite-field calculations. The refining parameter of our 
algorithms is the so-called bond dimension which con¬ 
trols the amount of entanglement in the tensor network, 
and which we consider up to = 15. In algorithm 2, Our 
PEPS is built from a unit cell of 9 sites together with 6 
simplex tensors. Translation symmetry may be preserved 
as in the 3-site pattern of Eig. 1(a) (point group symme¬ 
try only is broken), or spontaneously broken as in the 
\/3 X a/3 superstructure, i.e. 9-site pattern of Eig. 1(b). 
The characterization of the phases is possible by checking 
the local longitudinal and perpendicular magnetizations 
as a function of the magnetic field, as well as the local 
energies (Si-Sj) between all nearest sites within the unit¬ 
cell. To estimate expectation values of local observables, 
we consider here an approximate mean-field-like effective 
environment around the unit cell. 

Phase diagrams.- As we apply an external magnetic 
field, the system goes through several field-induced 
phases summarized in Table I. In Eigs. 2(a - d) we show 
the computed phase diagrams, presenting the (reduced) 
longitudinal magnetization as a function of the field for 
S' = ^, 1, I and 2. In the plots we can see several magne¬ 
tization plateaus, corresponding to incompressible phases 
with different symmetry breakings, and where the mag- 
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(a)S=1/2 (b)S=1 



FIG. 2: [Color online] Longitudinal magnetization (normalized w.r.t. its value at saturation) versus the external magnetic 
field for (a) S = (b) S' = 1, (c) S' = | and (d) S = 2, with bond dimension D = 10. The |-plateau (in red) is present for 
all spin values, and within our accuracies it is compatible with either a nematic or a solid phase. The dark blue phase, with a 
commensurable value of the magnetization (1 — ^), corresponds to the 2-magnons phase. The green part corresponds to the 
1-magnon phase with magnetization (1 — ^). Both the 1- and 2-magnon states are solid phases. For S = |, an additional 
incompressible phase (in light blue) exists at | of the total magnetization, which is also a solid phase but with coexisting VBC 
order (see text). Moreover, for integer spin values a spin gaped phase is identified in brown for low value of the field, with zero 
magnetization and corresponding to a simplex solid phase for S' = 1 and, to the best of our accuracy, to a topologically trivial 
spin liquid for iS = 2. Yellow regions are compressible supernematic or supersolid phases. Note that a small jump, highlighted 
by a bond dimension T) = 15, in the S' = 1 magnetization curve (at low field) signals a first order transition between two such 
compressible phases (see Ref. 9). Grey regions could not be fully characterized within our approach. The phases are labelled 
according to the notations in Table I. 


netization attains a commensurate value. Three of these 
incompressible phases are common to all the local spin 
values S. These are the plateaus at |th, (1 — ^)th and 
(1 — ^)th of the saturated value of the magnetization. 
In what follows we discuss some of the observed features 
in these diagrams. 

(i) ^-plateau: all the studied values of the local spin 
show a plateau in the magnetization at one third of the 
saturation value. A full characterization of this phase 
is subtle. Two competing phases, the Nematic and the 


Solid 1 phase, are degenerate within the accuracy of our 
approximations. While both of them are 3-fold degen¬ 
erate and two sites are always equivalent on every tri¬ 
angle, the Solid 1 phase breaks translation symmetry 
(see Fig. 1(c)), whereas the Nematic phase breaks the 
Cq point group symmetry. 

(ii) 1-magnon (1 — ^)-plateau: this plateau also ap¬ 
pears for all values of S. It can be understood exactly by 
a one-magnon picture (see Ref. [19]), and has the symme¬ 
tries of the Solid 1 phase (see Fig. 1(c)). In the x 
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superstructure, each independent magnon (flipped spin in 
a polarized environment) is localized on an hexagon with 
Cq symmetry. The one-magnon state is, in fact, also 
an exact PEPS with bond dimensions D = 6. The two 
boundaries of the phase are first order phase transitions, 
ending in a direct jump to the saturation value. 

(hi) 2-magnon (1 — ^)-plateau: as before, this plateau 
is also common to all values of the spin S. It was first un¬ 
derstood in terms of a \/3 x \/3 superstructure of (quasi¬ 
independent) two-magnon states [19] (Solid 1). Our nu¬ 
merical optimization, however, shows that the Cq sym¬ 
metry is lowered down to C2, while the two reflexions 
au and ay are preserved. Such a phase is called Solid 2, 
according to Table I (see Fig. 1(d)). Although the Solid 
1 and Solid 2 phases might be quasi-degenerate, the two- 
magnon phase with the Solid 1 pattern was not stabilized 
within our optimization scheme. 

(iv) VBC plateau: independently of the common in¬ 

compressible phases of semi-classical origin described 
above, we also And a specific “quantum” plateau for 
S' = ^ at niz = This phase corresponds to a 

VBC-Solid (VBCS), where all symmetries are broken 
except for the (7^^-reflexion, thus being 18-fold degen¬ 
erate (see Fig. 1(e)). The corresponding pattern can 
be schematically viewed as polarized spin-1/2 forming 
aV3 X a/ 3 superstructure (solid order similar to the S2 
phase) superposed with (phase-locked) dimerized bow-tie 
chains running in the au direction (VBC order). Note a 
grand canonical DMRG calculation [6] also revealed this 
plateau but associated it to a topological Z3 spin liquid 
with no symmetry breaking, in contrast to our findings. 

(v) Other features: furthermore, for integer spin values 
we identify a spin gapped phase at low enough magnetic 
held. This phase has zero magnetization and corresponds 
to a simplex solid phase for S' = 1, whereas for S = 2 
it corresponds, to the best of our accuracy, to a spin liq¬ 
uid with no symmetry breaking. We believe the simplest 
scenario allowed by the Lieb-Schultz-Mattis-Hasting the¬ 
orem [20] occurs: the S = 2 spin liquid is topologically 
trivial, its GS being non degenerate. In that case, it 
belongs to the same class as (and can be adiabatically 
connected to) the S = 2 Affleck-Kennedy-Lieb-Tasaki 
(AKLT) state [21], a universal resource for measurement- 
based quantum computation [22], for which we also found 
a spin gap [9]. However, this conclusion disagrees with 
the coupled cluster calculation [10] which states that the 
GS of the spin S' = 2 of the zero-field Heisenberg model 
has a non-zero magnetization. As shown in Fig. 2, we also 
And, for all the spin values, compressible supernematic or 
supersolid phases which bear the same spontaneous lat¬ 
tice symmetry breaking as their nearby solid phases. 

Quantum-elassieal erossover.- Large values of the spin 
S tend to quickly suppress quantum fluctuations in the 
system. In order to study this effect, we have analyzed 
the rate at which the |-magnetization plateau disappears 
as S gets larger. Our results are for bond dimension 
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FIG. 3: [Color online] Properties of the |-plateau phase, (a) 
Normalized critical magnetic fields Hci (circles) and Hc 2 
(squares) versus the inverse of the spin. Values for bond di¬ 
mensions D = 10 (blue) and D = 5 (red) are shown, as well 
as results from exact diagonalization and from perturbation 
theory for classical spins [23] . (b) Magnitudes of the longitu¬ 
dinal local magnetizations along (mJ) and opposite {m ~) to 
the field, computed with T) = 10. Dashed lines are polynomial 
fits in 


up to T) = 10 and spin up to S' = 2, and these are 
presented in Fig. 3. We plot the two critical lines Hci 
and Hc 2 in Fig. 3(a). For completeness we also show 
the results of an exact diagonalization for small clusters 
from Ref. [23]. For S = oc, i.e., classical spins, we expect 
the magnetization curve to become smooth and, there¬ 
fore, both critical lines end up at the same classical value 
Hci/S = HC 2 IS = 2. We indeed observe that, moving 
away from S = |, the size of the plateau shrinks quickly 
with increasing ^. We also plot the local magnetizations 
for the I-plateau versus for = 10. As shown in 
Fig. 1(b), there are two non-equivalent on-site magneti¬ 
zations, m/ > 0 (in the held direction) and m~ < 0 (op¬ 
posite to the held direction), such that = |- 

In the classical limit, the quantum fluctuations disappear 
and mt = —m~ = 1. 


Conelusions and outlook.- Here we have studied the 
spin-S' KHAF in the presence of a magnetic held up to 
S = 2 using tensor network methods based on an iPEPS 
algorithm and an enlarged supercell. We have charac¬ 
terized a host of (zero temperature) quantum phases, 
some of semi-classical origin, some purely quantum. Our 
results confirm and extend in several directions previ¬ 
ous results in the literature : we show that the well- 
known S = 1/2 ruz = I and ruz = 1 — ^ plateaus 
extend to larger spin S. The spin-1 rUz = 0 sim¬ 
plex solid plateau remains stable within our new en¬ 
hanced variational space. However, we found that the 
spin-1/2 rUz = 1/9 plateau is stabilized by spontaneous 
breaking of translation symmetry, in contrast to previous 
claims. Also we found a less-symmetric pattern for the 
rUz = 1 — ^ two-magnon plateau than proposed previ¬ 
ously. One important finding of the present work is also 
the first full characterization of all the compressible (i.e. 
superfluid) phases between the plateaus. 
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Under completion of this work, we became aware of 
a related work showing similar findings on the Husimi 
lattice [24]. 
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I. OPTIMIZATION METHOD: SIMPLE UPDATE 
A. Description 

Tensor Networks in 2 dimensions in the thermodynamic limit, called infinite Projected Entangled Pair State 
(iPEPS), can be optimized using Eull Update or Simple Update schemes. The Eull Update is more accurate but 
requires the calculation of the surrounding system of the unit cell - called environment - at every step. The calcula¬ 
tion of the environment can be performed using a Matrix Product State (MPS)-based approach [3], a Corner Trans¬ 
fer Matrix Renormalization Group (CTMRG) [3] or a Coarse Graining Tensor Renormalization Group (GGTRG) 
method [4] (or any related technic). In contrast, in the Simple Update, the tensor optimization used in this study, 
the computation of the environment is not done during the optimization stage. The resulting optimized tensor may 
have a lower accuracy, but the method allows to go to higher values of the bond dimension D. 

In this paper we use two types of updates. The first one (’’algorithm 1”) is based on the mapping of the kagome 
lattice to the square lattice shown in Eig. 1.3, and is the one used in the calculation of the energy at zero field for 
the spin-1/2 case (see below). In the second one (’’algorithm 2”), we use the imbedded cluster of nine sites depicted 
in Eig. LI, which is reliable if quantum fluctuations do not extend much beyond the unit cell (it becomes exact if 
one considers a Bethe lattice). This iPEPS representation considers a large cluster (the focused triangle with its nine 
neighboring triangles as shown in Eig. 1.2), which has proved convenient for the calculations at finite field. Given its 
efficiency, in this paper we use this type of update to compute the phase diagrams. 





FIG. I.l: Tensor Network of 9 sites per unit cell, composed by 9 site tensors of dimension dD^ and 6 simplex tensors, located 
on the triangle, of dimension . 






FIG. 1.2: Cluster used in the Simple Update approximation. The red, green and blue site tensors picture the A, B and C 
tensors. The simplex tensors with circle (square) ends depict the down (up) triangle simplex tensors. 


B. Algorithm 1 

The mapping in Fig. 1.3 preserves the locality of interactions in the model (one has, at most, nearest-neighbor 
interactions also in the square lattice). The system is thus amenable for simulations with a standard iPEPS algorithm 
on a square lattice, with a squared physical dimension. If a 2-site unit cell is used, then the calculations are, essentially, 
an improvement upon a 6-site unit cell in the original kagome lattice. The improvement comes from the fact that the 
three sites that are coarse-grained are treated by a single tensor instead of three. In our case, here we used also the 
simple update scheme. 


C. Algorithm 2 

The optimization of the tensors in the Simple Update is based on A matrices which are located on the external 
bonds of the triangles. In one dimensional system, these A matrices contain the singular values of a decomposition 
of two neighboring sites, and connect these two sites. In two dimensions, they are an approximation of the outside 
system, whereas the simplex tensor connect the three sites inside the focused triangle. 

In the following, the algorithm [1] for one optimization of the tensors on down triangles around simplex tensor 
RV[l] 

is described (see Fig. 1.4). 

1. First, we contract the T tensor on the triangle, as shown in Fig. 1.4(a), 

rj.SA,SB,Sc _ MSa-d[1]SbMI]Sc ),AA[1]. AB[1]. AC[1] 

l,m,n 


2. Then, we apply the imaginary time evolution operator U = exp on T: 
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FIG. 1.3: Mapping scheme from the kagome to the square lattice. We combine the three orange marked spins together and 
end up with a square lattice. One single link in the square lattice corresponds to two links in the kagome lattice. Since all 
nearest neighbours in the kagome lattice are mapped to either one single site or nearest neighbours in the square lattice, the 
interaction remains local. 



FIG. 1.4: T tensor contraction, including diagonal matrices on the down triangle. 
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where S' is a vector, whose components are Sa, Sb and Sc, and Sr 1 (Fig. L4(b)). 

3. Now, we have to decompose the V tensor into the new tensors and This step is called a 

Higher-Order Singular Value Decomposition (HOSVD). We give the example for the tensor. We diagonalize 
the dD X dD matrix W = as shown in Fig. L4(c). 




(5a, i) iO 


SB,Scd,k 


Sa,Sb,Sc 


4. We keep the D largest eigenvalues of the decomposition in and is the corresponding dD x D matrix. 
The new tensors are -yj/ (Fig. F4(d)) and (where > 0 because = W). 

After the updates on the three sites, we can compute the new simplex tensor as 

Klk = E E (bT.O;0* (bl.m).)* (b^cn);.)* 

{Z,m,n}{5A,5B,5c} 


The renormalization on up triangles is identical, by doing the change A taking 

and Then, we perform two cyclic permutations on the labels [1 ^ 2 ^ 3] in order to optimize all the nine 

site tensors and the six simplex tensors. In order to reach the ground state, one can use a succession of several time 
steps from Sr = 10“^, 10“^, 10“^, 10“^, 10“^. For each time step, one can choose a accuracy criterion for which the 
evolution is stopped. In addition, one can stop the evolution if the number of time step is too large. Here, we used 
the normalization of the lambda matrices with a relative accuracy A = 10“^^. When the number of time steps was 
too large, we stopped the imaginary time evolution at 4000 iterations. In those cases, the order of the accuracy was 
10-^ 


II. ENERGIES 

Here we show some results for the ground state energy for spin 1/2 and zero field using algorithm 2 above. Our 
calculations in this case are for a PEPS with a 6-site unit cell in the kagome lattice, where effective environments are 
computed by Corner Transfer Matrix (CTM) methods, up to PEPS bond dimension D = 10. The results are shown 
in Eig.H.l, where the points are converged in the environment bond dimension. The dashed line corresponds to an 
algebraic fit, leading to asymptotic ground state energy Eq{D -d oo) ^ —0.436. Our numbers obtained using this 
approach are, in fact, very close to the ones obtained using the PESS approach, which is = —0.4364 for bond 

dimension D = 13 and a 9-site unit cell [1]. 

At finite field, we can show that, using algorithm 1 above, the expectation value calculation using only the lambda 
matrices (called simple calculation) is a good approximation compared to the full calculation^ which corresponds to 
the calculation of the full surrounding environment. We show this comparison in the case of the spin S' = 1 for the 
nematic {h = 2.6) and supernematic {h = 4.29) phases (Eig. H.2(a-b)) for which the unit cell is made of only three 
different sites. The environment was calculated using the MPS-based approach, with the couple {D^Dc}, where Dc is 
the environment bond dimension, {3,18}, {4,23} and {5,30} for the nematic phase and {3,29}, {4,29}, {5,31} and 
{6,30} for the supernematic phase. 


III. REDUCED MAGNETIZATION 

As we apply an external magnetic field, the systems go through several phases. In order to identify these different 
phases, we calculcate local expectation values using the procedure describes in the Supplemental Material [2]. Eirst, we 
compute the reduced longitudinal magnetization, as presented before, on each site = {S^)/S, where the label fa G 
[1, 9] runs over the nine sites. The average of the magnetization along the magnetic field m^/9 is especially 

useful for identifying the several magnetization plateaus. The local reduced magnetization s along the m agnetic 
field together with the local reduced magnetizations transverse to the magnetic field /S give 

information about the symmetries of the phases. The transverse magnetization is shown in Eig. HI.l. The SU{2)- 
symmetry breaking is characterized by a non-zero local magnetization in the transverse plane. Since the plateau 
phases preserve this spin symmetry, all the local transverse magnetization are zero, as shown in Eig. HI.l, for every 
spins. Moreover, in every phases, but the gray regions, the global SU (2)-symmetry is preserved, which is characterized 
by = 0 
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FIG. ILl: Energy of the spin-1/2 case, h=0, for a PEPS with a 6-site unit cell in the kagome lattice, versus the bond dimension 
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EIG. IL2: Gomparison between the simple and the full calculation of the spin-1 case for the nematic (a) and the supernematic 
phase (b). 


In addition, we calculate the eigthteen local energies Cij = (Si • Sj) for i and j nearest neighbours. This density is 
pictured in the Fig. l(c-e) of the main paper and is necessary to characterize the Valence Bond Crystal phase. 
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FIG. IILl: [Color online] Reduced transverse magnetizations (normalized w.r.t. its value at saturation) versus the external 
magnetic field for (a) S' = (b) S' = 1, (c) S' = | and (d) S' = 2, with bond dimension D = 10. The plateau phases 

(in red, light blue, dark blue and green) and the gaped phases for spin S = 1 and S = 2 are characterized by a zero local 
transverse magnetization for every site. Yellow regions are compressible supernematic or supersolid phases. Although the local 
magnetizations are non-zero, the total transverse magnetization is zero = 0. Grey regions could not be fully 

characterized within our approach as show the transverse magnetization behavior. The phases are labelled according to the 
notations in Table I and the colors correspond to the Fig.2. 
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